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Abstract 

The crossover to Rouse- like behavior for the self-diffusion constant D, the 
viscosity rj, and the equilibrium structural statistics of n-alkanes (6 < n < 66) 
is studied numerically. For small n the chains are non-Gaussian and the mean 
squared end-to-end distance (R 2 ) is greater than 6(Rq), where (Rq) is the 
mean squared radius of gyration. As n increases, (R 2 }/(Rq) — ► 6(1 + a/n), 
where a depends on the interaction model. At constant density, the Rouse 
model is used to extract the monomeric friction coefficient C and the viscosity 
r] independently from the diffusion constant D and the longest relaxation 
time tr. Cd extracted from D is nearly independent of chain length while £ T 
obtained from tr is much larger than £d for small n. The viscosity measured 
in a non-equilibrium molecular dynamics simulation is closely approximated 
by the value of r\ determined from tr while r\ inferred from D is smaller for 
small n. For n > 60, the two estimates for both £ and r] agree as predicted 
from the Rouse model. D calculated from three interaction models is studied 
for increasing n and compared to experimental data. 
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I. INTRODUCTION 

It is well known that dynamics of a melt of unentangled polymer chains satisfies the 
Rouse model [|1] exceedingly well. For chains which are shorter than an entanglement length 
n e , the self-diffusion constant and viscosity for a melt of chains of length n are observed 
experimentally to scale as n~ l and n, respectively, as predicted by the Rouse model. 
Computer simulations on simple bead spring models [§J have shown that not only does 
the Rouse model predict the correct scaling dependence on chain length n, it also gives 
quantitatively correct results. Why this simple model, which was first proposed |l| to model 
the dynamics of dilute polymer solutions, works so well for a polymer melt remains somewhat 
of a mystery. 

The crossover from Rouse-like behavior for unentangled short chains to reptation-like 
motion for entangled chains has received considerable attention in recent years, while the 
crossover from the behavior of very small molecules to the Rouse regime has received rela- 
tively little attention. Understanding the latter crossover can provide information on sub- 
tleties of the Rouse model and at the same time increase knowledge on the dynamics of a 
technologically important class of molecules, normal alkanes. There have been a number 
of studies of the dynamics of short alkanes, but to our knowledge there have been no sys- 
tematic studies which examine the dependence of either the diffusion or viscosity on chain 
length. One can also examine the crossover of the mean squared radius of gyration (Rq) 
and end-to-end distance (R 2 ) as a function of n. For Gaussian chains, (Rq) and (R 2 ) satisfy 

(R 2 ) = 6{R%) = nb 2 (1) 

where b is the effective bond length. Consider the simple case of normal alkanes, though the 
same is true for any flexible chain molecule: for small n, the chains are clearly non-Gaussian 
and (R 2 ) > 6(R G ). For normal n-alkanes, crossover to Gaussian chain statistics occurs for 
n greater than 100 |||6],|7|] . Since the Rouse model is based on the fact that the chains are 
Gaussian, one would expect that crossover length is a minimum chain length for the Rouse 
model to hold, though this has not previously been checked systematically. 
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In the Rouse model, the excluded volume interactions and the hydrodynamic interactions 
are disregarded and the polymer is treated as a collection of beads connected with a harmonic 
spring with spring constant k = 3k B T '/b 2 '. Here T is the temperature and k B is Boltzmann's 
constant. The dynamics of the chain is then modeled with a Langevin equation with a 
monomeric friction coefficient (. The resulting set of equations can be solved by transforming 
to normal coordinates. The resulting self-diffusion constant is |§ 

k B T , \ 

D = -V- 2 

nQ 

The rotational relaxation time tr of a polymer is usually defined by the longest relaxation 
time of the end-to-end autocorrelation function (R(i) • R(0)) and expressed as f8|, 

TR ~ 3^T- (3) 
The effective bond length b can be eliminated in favor of either (R 2 ) or (Rq) using eq. [I], 

C^ 2 ) 2(n(R 2 G ) 

R 37T 2 k B T 7T 2 k B T ' 1 ' 

The viscosity rj of a melt can be determined from the relaxation modulus G(t), which can 
be written as a sum over the Rouse modes of a chain. In terms of r R ||, 

TT 2 pRTr R 

^ = ^2iW-' (5) 

where M is the molecular weight, p is the density, and R is the gas constant. This gives 
the well known result that for Rouse chains, the viscosity scales linearly with n for constant 
monomeric friction coefficient. This simple scaling does not hold at constant temperature 
and pressure, since p and therefore ( increase with increasing n. This results in a stronger 
dependence on n, namely 77 ~ D^ 1 ~ n^ 1 ' 8 for unentangled polyethylene chains for T = 
448 K §. 

To study the crossover from small molecule behavior to Rouse-like chains, it is convenient 
to rewrite the monomeric friction coefficient ( and viscosity rj in several different forms. From 
eq. for the diffusion constant D, we can write 



CD = ~nD' (6) 

From eq. ^ for the rotational relaxation time tr, we have two measures of ( depending on 
whether we use (R 2 ) or (Rq) to eliminate b, 

37T 2 k B Tr R 



CAR) 



(ARg) = -rrm 



n(R 2 ) 

(7) 

n 2 k B Tr R 
2n(R 2 G ) ■ 



The viscosity can also be written three different ways. The first, given by eq. [5], we refer to 
as r\ T . Two additional expressions for rj in terms of D, can be obtained by substituting each 
form of eq. [| for tr in eq. |S| and eliminating ( in favor of D using eq. |2|, 

P RT(R 2 ) 
Vd{R) ~ -36MD- 

(8) 

These three expressions for both ( and rj are equivalent for Rouse chains, so that any 
differences are a clear signature of deviations from Rouse behavior. 

To investigate the crossover to Rouse-like chains with increasing chain length, we car- 
ried out a series of equilibrium (EMD) and non-equilibrium molecular dynamics (NEMD) 
simulations for normal alkanes with 6 < n < 66 at the density for high molecular weight 
polyethylene. Constant density was studied due to the large corrections to Rouse behavior 
observed experimentally for unentangled chains below the entanglement length n e at atmo- 
spheric pressure. The EMD simulations were used to calculate both ( and rj as predicted 
by the three expressions for each discussed above. The results for rj, as predicted by the 
EMD calculations, were compared to the viscosity calculated directly from the NEMD sim- 
ulations. We also carried out EMD simulations on alkanes in this size range employing three 
different interaction models discussed recently in the literature. These were performed at 
experimental density for atmospheric pressure and the diffusion constant D as a function 



of chain length was compared to experiment. All sets of EMD simulations were used to 
study the behavior of the equilibrium structural statistics (Rq) and (R 2 ) with increasing 
chain length. To efficiently extend this comparison out to sizes larger than n = 66, we also 
carried out simulations using the rotational isomeric state (RIS) model for chains in the 
range 34 < n < 2000. 

In Sec. II we summarize the simulation models and the methodology used. The equi- 
librium results for the static structural properties and their behavior with increasing n are 
presented in Sec. III. Results for the dynamics are presented in Sec. IV where we discuss the 
crossover for the monomeric friction coefficient and viscosity as a function of chain length 
n and compare diffusion results to experiment and predictions based on the Rouse model. 
We briefly summarize our main conclusions in Sec. V. 
II. SIMULATION MODEL AND METHOLOGY 

In order to study a wide range of chain lengths, an united atom (UA) model was used to 
simulate the linear alkanes. In this model, the hydrogen atoms are not explicitly simulated. 
Rather, they are grouped with the carbon to which they are bonded into a single particle 
or a united atom. We do, however, distinguish between different CH n groups. These groups 
or united atoms interact through bonded and non-bonded forces. The bonded interactions 
are represented by constraint forces which keep intramolecular nearest-neighbors at a fixed 
distance, a bending term 

W) = - e b f ( 9 ) 

where 9b is the equilibrium angle between successive bonds, and a torsional term 

\4(0)=5>cos*(0) (10) 

i 

characterizing preferred orientations and rotational barriers around all non-terminal bonds. 
The non-bonded forces are described by Lennard- Jones (LJ) interaction sites located at the 
position of each carbon atom center of mass. The LJ potential is defined by 

V LJ (r) = Ae\C-f 2 - i^f] (11) 
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This non-bonded interaction is between both intermolecular and intramolecular sites except 
for those which are separated by less than four bonds and which therefore interact through 
one or more of the bonded interaction terms. The LJ interaction is usually truncated at a 
distance r c and the potential shifted so that Vlj{t c ) = 0. 

In our earlier studies of the diffusion [II], DJ and viscosity |T2j of linear and branched 



alkanes, we studied two models for the interactions. The first was the symmetric UA model of 
Siepmann et al. fl3| , |l4|1 which was developed to describe the vapor-liquid coexistence curves 
of the n-alkane phase diagram ||15|| . We refer to this as model A. The second, which we 



denote as model B, was the asymmetric united atom (AUA) model of Padilla and Toxvaerd 
|16| , p!7[| which was optimized to describe the static and dynamic behavior of short n-alkanes 
(5 < n < 10) at moderately high temperatures and pressures. In this latter model, there is a 
displacement between the centers of force of non-bonded interactions and the centers of mass 
of the united atoms, hence the name. Paul et al. |18[ recently introduced a new, optimized 



symmetric UA model which they claimed was better able to reproduce the properties of 
melts of linear normal alkanes. The primary differences between their model and that of 
Siepmann et al. JT3"| , |il| was the use of a slightly larger a and a different torsional potential; 
we refer to our implementation of this last parameterization as model C [ T9 ]. 



The intramolecular interaction parameters for the three models are listed in Table |. 



The parameters for the torsional potential of model A are from Jorgensen et al. [20|. Paul 



et al. Jl8| presented results for two torsional potentials. The first one they studied gave 
a diffusion constant which was 45% larger than the experimental result for n — C^H^ at 
T = 400 K. They subsequently modified the torsional potential V t ((f)) to increase the trans- 
gauche barrier from 3.0 to 3.3 kcal/mol. Since this latter model gave a better estimate of the 
D, though still 30% too large for n — C i4: H 9Q , we used their modified version of V t ((j)) here. 
For comparison, the torsional potentials for the three models are shown in Fig. [I]. Note that 
the gauche energy is about the same for models B and C (0.5 kcal/mol) and about half that 
of model A (0.9 kcal/mol). This difference will influence the thermodynamic flexibility of 
the molecules and the persistence length for long chains. Models B and C are expected to 



result in more thermodynamically flexible molecules than model A. The dynamic flexibility 
is related to the size of the trans-gauche barrier, which is larger for model B than for models 
A and C. Thus, models A and C lead to more dynamically flexible molecules than model B. 

The Lennard- Jones parameters are listed in Table |I| for all three models. Note that the 
interaction strength e is the same for models A and C. The effective diameter for the united 
atoms is different for the three models, with model C having a larger effective diameter than 
model A. For model B, a direct comparison based on the value of a is not appropriate due 
to the asymmetric nature of united atoms. While the a for model B is the smallest, the 
displacement between the center of force of non-bonded interactions and the center of mass 
of the united atom, results in an effectively larger united atom than represented by a. The 
interaction parameter between a CH 3 and CH 2 site is given by the Lorenz-Bertholot rule, 
v /e CH 3 ecH 2 - We used a 10 A cutoff for the (shifted) LJ potential. 

All of the MD results presented here were obtained from constant volume, constant 



temperature simulations. The velocity rescaling algorithm of Berendsen et al. Ell] was used 



to control the temperature in the EMD simulations while the Nose-Hoover [E^] algorithm was 



used for the NEMD simulations. The equations of motion were integrated using the velocity 



Verlet algorithm P2[ with a 5 fs time step for both the EMD and NEMD simulations. Bond 



lengths were kept constant using the RATTLE algorithm ^2[^3j. Equilibrium structural 



statistics, the self-diffusion constant D, and rotational diffusion time tr were determined 
from EMD simulations while the viscosity rj was determined from NEMD simulations. The 
NEMD results reported here were obtained in the Newtonian regime, for a shear rate 7 = 
8.70 x 10 -4 ps _1 , which satisfies 7 < r^ 1 for all n studied. 

One set of EMD simulations was done at constant density for a series of chain lengths 
n. The density used was the bulk density for high molecular weight polyethylene, p = 
0.766 g/cm 3 and the simulations were carried out at (T = 448 K). Selected chain lengths 
at this state point were also studied with NEMD simulations. For small n, this state point 
corresponds to a high pressure state. In all of the constant density simulations (EMD 
and NEMD), interaction model A alone was employed. Another set of EMD simulations 



were performed using the experimental densities at atmospheric pressure (P ~ 0.1 MPa) 
for four chain lengths n. For each of these chain lengths, simulations were run using all 
three interaction models A, B, and C. These experimental density simulations were run 
at T = 400 K in order to remain far below the vapor-liquid critical temperature T c for 
small n, which for example is 507 K for n-hexane and 617.5 K for n-decane P^] . Since we 
were interested in studying the effect of chain length on diffusion and viscosity, we chose 
temperatures which were higher than the melting temperature T m for large n. For linear 
polyethylene, T m ~ 400 K. The lengths of the runs and the number of molecules used for 



each case are presented in Table |T| and Table [TV. 

To examine equilibrium structural statistics for n > 66, we used a continuum version of 
the RIS model for chains up to length n = 2000. There have been several studies |5|,|6|,[7|, |25|j26 



which have shown that by proper choice of the cutoff for the non-bonded interactions, the 
radius of gyration and end-to-end distance of an isolated chain compare very favorably 
with that of the melt. Here we followed the standard RIS model for polyethylene, which 
includes all bonded interactions and the 1-5 intramolecular Lennard- Jones potential (these 
are the "second-order interactions" of RIS theory and are known to be essential in order 
to reproduce the experimental characteristic ratio of the polyethylene chain). The main 
difference between the present model and the standard RIS model is that in the present 
implementation the molecule can explore the entire torsional potential surface, rather than 
being confined to a set of its minima. To allow comparison with the constant density EMD 
simulations, the simulations were done at T = 448 K. For short chains (n < 100), one can 
use Langevin dynamics simulations in which each isolated chain is coupled weakly to a heat 



bath. However for large n, the pivot algorithm [[2t[ is considerably more efficient. Most of 
the results for the RIS simulation were obtained using the pivot algorithm, though for small 
n both methods gave comparable results as expected. 
III. EQUILIBRIUM PROPERTIES 

The structural properties for the molecules as expressed by the mean squared radius of 
gyration (Rq), the two largest eigenvalues of the mass tensor I] and Z|, and the mean squared 
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end-to-end distance (R 2 ) are collected in Tables [TTT] and |V| for all the simulations. Results 



for (R 2 )/Q(R 2 G ) are shown in Fig. @a for model A at T = 448 K. Data for both the EMD 
(p = 0.766 g/cm 3 ) and the RIS simulations are presented. Note that for small n < 16, this 
ratio increases with increasing n. For small n, the chains are relatively rigid. As such, the 
change in chain dimensions with increasing n is similar to what can be calculated for a linear 
string of uniform beads, where {R 2 ) /6{Rq) goes from ~ 1 for very small n to 2 for infinite 
n. Above n > 16, flexibility of the chain becomes significant and (R 2 ) /6(Rq) decreases with 
increasing n. A similar result was found by Baschnagel et al. f| and Brown et al. ||. The 
EMD results show that, for the values of n studied, the chains are clearly extended and 
non-Gaussian as (R 2 ) > 6(Rq). For this model and temperature, the crossover to Gaussian 
statistics occurs for n much greater than 100 in agreement with Baschnagel et al. ||, Brown 
et al. || and Paul et al. [^]. Figure |a shows that the values of (R 2 ) /6{Rq) predicted by 
RIS and EMD simulations agree within our statistical uncertainty. For the longest chain 
length studied, n = 2000, (R 2 ) /6(Rq) is roughly 1% greater than the prediction for Gaussian 
chains for model A. A second order curve has been constructed on Fig. 0a representing the 
best fit of all the data for the RIS simulations and the data from the EMD simulations at 
n = 66,44,36,30,24, and 16. For model A, ((R 2 )/6(R 2 G ) = 1 + 11.8/n - 93.4/n 2 ). From 
this curve one can form an approximation of the chain length needed to reach Gaussian 
statistics within a given degree of error. For instance, the deviation from Gaussian statistics 
is less than 2% for chain lengths n > 600 for model A. This is considerably larger than 
for a flexible chain consisting of beads connected with a spring jp but with no torsional or 
bending forces, in which (R 2 ) /6(Rq) = 1 for n > 10. 

From our previous studies for small n (n < 24), we expect that the actual crossover to 
Gaussian chain structural statistics will depend somewhat on the specific model, particularly 
through the parameters of the torsional potential. As discussed in Sec. II, the torsional 
potential for model A has a larger gauche bond energy than that used in models B and 
C (about 0.4 kcal/mol larger). This generally means that models B and C result in more 
thermodynamically flexible, and therefore more compact, molecules than model A. Data in 
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Table [TV| confirms these expectations; furthermore, the ranking of (R 2 ) for the three models 



is constant with n: (R 2 ) A > (R 2 )c > (R 2 )b- Comparing data for model A and model C, 
the influence of the gauche bond energy on (R 2 ) is evident. This comparison is valid as 
models A and C represent the same implementation of the united atom model, only with 
different parameters. Furthermore, their largest difference in parameters, for this study of 
equilibrium properties, is found in their torsional potentials. Thus the effect on molecular 
size can be directly attributed to the difference in gauche bond energy. Comparing (R 2 ) 
data for models B and C, we note that the difference in molecular size for the two models 
is larger than what one would expect due simply to the slight difference in the gauche bond 
energy. This is sensible since model B represents a fundamentally different parameterization 
of the united atom model than models A and C. This indicates that torsional potentials are 
not directly transferable between UA and AUA models. From observation of the data in 
Table [TV| , we anticipate that the crossover to Gaussian structural statistics with increasing 
n will occur first for model B, then for model C, and last for model A; Fig. |2|b gives evidence 
of this. This figure is the same as Fig. |2]a but with data from the EMD simulations using the 
three models at T = 400 K and experimental density at atmospheric pressure. While Fig. ^|b 
shows only the beginning of the crossover, the ranking predicted is evident at this stage. 
The RIS simulations for the three models confirm the order of the crossover to Gaussian 
statistics. 

Differences in the torsional potential also give rise to a difference in the expansion co- 
efficient k = d(ln(R G ))/dT between the models. In a previous publication [10], we found 



that for tetracosane around room temperature, k = -1.1 ± 0.2 and -0.6 ± 0.3 x 10~ 3 K" 1 
for model A and B, respectively. These results were for simulations run at the experimental 
density for P = 0.1 MPa for each temperature. This data was compared to the experimental 
value for bulk polyethylene k = —1.07 x 10~ 3 K _1 obtained from SANS measurements [|28| . 



We can also consider previous data for (Rq) of n-hexadecane |L2| obtained at T = 298, 323, 
and 373 K in conjunction with (Rq) obtained at T = 400 K presented in this paper. From 
these, we can calculate k = —0.8 ±0.2 and —0.6 ±0.1 x 10~ 3 K _1 for n-hexadecane in model 
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A and B, respectively. These results are not presented to further evaluate the two models so 
much as to demonstrate that there is a very slow convergence of molecular sizes for models A 
and B with increasing temperature. In this respect, note that the average radius of gyration 
of n-hexadecane using model A at T = 400 K is still larger than the value for model B at 
298 K flr0"|. While temperature can influence the position in n of the crossover to Gaussian 



statistics, these observations show that the model used has a much greater influence than 
temperature. 



Boothroyd et al. [gg] have reported on experimental data for the radius of gyration of 
a polyethylene chain in a melt. Specifically, they quote the radius of gyration for samples 
of molecular weight My/ = 32000 and 50000 in the temperature range of 380 K to 480 K. 
Dividing the experimental (R 2 G ) by M w for M w = 32000 and 50000 at 448 K, we obtain 0.23 
and 0.19, respectively. In Fig. |3], (Rq)/M w is plotted versus 100/M^ 2 . Data obtained from 
both the EMD and RIS simulations at T = 448 K are shown for model A. In addition, a 
linear fit to all the RIS data and the EMD data for 30 < n < 66 is also shown. Extrapolating 
this line to M\y — > oo gives a value of 0.28 for model A which is considerably larger than the 
experimental results for polyethylene. Similar fits to the RIS model using the parameters 
for models B and C give (Rq)/M w = 0.18 and 0.24, respectively, in better agreement with 
experiment. 

IV. DIFFUSION AND VISCOSITY RESULTS 

Results for the self-diffusion constant D, the rotational diffusion time tr and the shear 



viscosity 77 are collected in Tables [TTI| and [TV] for all simulations performed. The self-diffusion 



constant is obtained from the slope of the mean-squared displacement of the molecular 
center of mass, averaged over all the molecules (A mo i) in the system and all available initial 
times. In taking the average, the sum of the mean-squared displacements for the iV mo i is 
divided by A mo i — 1 to account for the fact that the center of mass of the system does 
not move. To calculate tr, we consider the first-order angular correlation of the longest 
principal axis of a molecule's ellipsoid of inertia [|Tt]]. At intermediate times, this correlation 
is well described by a simple exponential relaxation with a time constant tr as shown 
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in Fig. |] for n-tetracosane and n-hexahexacontane. Equivalently, one can measure the 
autocorrelation function (R(i) -R(0)) for the end-to-end distance []18| and fit this curve with 
the corresponding formula from the Rouse model || to obtain a measure of tr. Both of 
these correlations can be calculated from EMD. The viscosity was determined directly from 
NEMD simulations. From our previous study [|T^] comparing EMD and NEMD simulation 
methods for determining the viscosity of n-alkanes, the NEMD simulation technique requires 
about half the cpu time required by EMD methods for comparable accuracy. For further 



details see ref. [12 



The diffusion constants for the three models studied are compared to experimental data 
29| for T = 400 K in Fig. |5|. These simulations were run at the experimental densities listed 



in Table [Iy| for P = 0.1 MPa. Note that for all of the simulation results as well as the 
experiments, the product Dn decreases with increasing n. This decrease is due to the fact 
that the density and therefore the monomeric friction coefficient ( increases with increasing 
n, causing D at fixed T to decrease faster than predicted by the Rouse model. Assuming 
that the change in Dn is due to change in the monomeric friction coefficient a fit to the data 
gives ( ~ n~ a8 . One also sees from Fig. and Table [TV] that model A overestimates D by 
approximately 15 — 50%, while models B and C fit the available experimental data quite well. 
Similar results for models A and B were found in our earlier studies [fTo|JTT|JT^l for n-decane, 
n-hexadecane, and n-tetracosane as well as for several branched alkanes. What is new in the 
current results is the agreement with experiment found for model C which is a symmetric 
united atom model. A symmetric united atom model is somewhat less computationally 
demanding to implement than an asymmetric united atom model. For identical length 
simulations we find that for 100 molecules of n = 44, the symmetric UA model requires 
approximately 10% less cpu time than the AUA model. 

The three estimates of ( are shown in Fig. |6| for constant density and T = 448 K. Note 
that (d determined from the diffusion constant is only weakly dependent on n and even very 
small chain lengths give a reasonably good estimate of the friction coefficient of unentangled, 
Rouse-like chains. The change in (d from n = 10 to 66 is only 10%. The significant change 
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in ( with n at constant pressure as exhibited by the diffusion data discussed above all but 
disappears at constant density. This is in contrast with ( T which depends much more strongly 
on n for small n\ furthermore, £ T > C,d |P-8| . The difference between ( t (Rg) and ( T {R) for 
small n is another representation of the deviation of equilibrium structural statistics from 
Rouse-like behavior. For n > 60, all three methods of determining ( agree within our 
statistical uncertainties as predicted by the Rouse model. 

Another attribute of Fig. |^ warrants discussion. At small n, (d decreases with increasing 
n and appears to reach a minimum around n = 24. This behavior is caused by an end 
effect: at small n, there is a large concentration of CH3 united atoms whose Lennard- Jones 
e parameter in model A is over twice that for CH 2 united atoms. An increase in e is known 
to decrease diffusion |H| thereby effectively increasing (j> To verify that this end effect 
caused the increase in ( D , the state points for n = 6, 16, and 24 were simulated but with the 
e = 0.93 kcal/mol for all the united atoms (i.e. the value for CH2 in Model A). The mass of 
all united atom groups was that of CH 2 so the volumes were adjusted to achieve the stated 
density. The results for are shown on Fig. |^ as filled symbols. The change in value at 
n = 24 is very small (4%) demonstrating that the end effect is nearly undetectable for large 
n as expected. The overall variation in is now also smaller. From the three data points 
presented, it can be seen that the minimum in (d with n still exists, however the minimum 
has shifted to n = 16, coincident with the position of the maximum in (Rq) (Fig. ^|a). This 
suggests that the origin of the minimum of (d in Fig. |6] is related to the chain flexibility. 

Finally, the three estimates for the viscosity are compared to our results from NEMD 
simulations in Fig. [F[ As observed in our previous study |12| , r\ T gives a very good estimate 
(within about 20%) for all n, while r/o is significantly smaller for small n. Since tr requires 
approximately ten times less cpu time to determine than rj, we believe that at least for 
unentangled chains, tr is an efficient way to predict 7] without having to measure 7] directly. 
For n 60, all the estimates of 77 and the NEMD calculation agree within our statistical 
uncertainty. 
V. CONCLUSIONS 
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The crossover to Rouse-like behavior for n-alkanes (6 < n < 66) as demonstrated by 
their self-diffusion constant D, viscosity i], and equilibrium structural statistics was studied 
numerically. For small n the chains were non-Gaussian and the mean squared end-to-end 
distance (R 2 ) was greater than Q(R 2 G ), where (Rq) is the mean squared radius of gyration. 
As n increases, (R 2 ) / (Rq) — > 6 but only for n significantly larger than 200. At constant 
density, the Rouse model was used to extract the monomeric friction coefficient ( and the 
viscosity i] independently from the diffusion constant D and the longest relaxation time tr. 
(d extracted from D was nearly independent of chain length while ( T obtained from tr was 
much larger than (o for small n. The viscosity measured in a non-equilibrium molecular 
dynamics simulation was closely approximated by the value of 77 determined from t r while 
the 7] inferred from D is smaller for small n. For n > 60, the two estimates for both ( 
and rj agree as predicted from the Rouse model. Diffusion at normal pressure as calculated 
from three interaction models was studied for increasing n and the results were compared to 
experimental data. Very good agreement with experiment was found for two of the models 
while the third consistently overestimated D. For all models and experimental data, D was 
found to depend more strongly on n than predicted by Rouse theory. 
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FIG. 1. Comparison of torsional potential around X — CH2 — CH2 — Y type bonds for models 
A, B, and C. 



FIG. 2. Mean squared end-to-end distance (R ) divided by the radius of gyration (Rq) versus 
n. For a Gaussian chain (R 2 ) /6(Rq) = 1. (a) T = 448 K and p = 0.766 g/cm 3 . Data from model 
A for bulk and RIS simulations, (b) T = 400 K and experimental density at P ~ 0.1 MPa. Data 
from models A, B, and C for bulk simulations. 



FIG. 3. Radius of gyration (R 2 G ) divided by mass M versus 100/M 1 / 2 for T = 448 K and 
p = 0.766 g/cm 3 , model A bulk and RIS simulations for models A, B, and C. 

FIG. 4. Time autocorrelation function of the orientation of the longest principal axis of inertia 
for n-tetracosane (n = 24) and n-hexahexacontane (n = 66) for T = 448 K. The solid line is for 
(n = 24) and the dashed is for (n = 66). 

FIG. 5. Diffusion constant D from models A (□), B (A), and C (o) for T = 400 K at the 



experimental densities listed in Table |Vj for P = 0.1 MPa. The experimental results (•) are from 
ref. @. 

FIG. 6. Friction coefficients £ versus n for T = 448 K. (o) is extracted from the self-diffusion 
constant, eq. |, while Cr(R) (□) and Ct(Rg) (A) are determined from the longest relaxation time 
r H) e Q- @- The filled symbol (•) is C,d obtained from simulations for chains consisting entirely of 
CH2 united atoms. 

FIG. 7. Viscosity i] T (o) determined from the longest relaxation time tr, eq. |j| and rfn(R) (Q) 
and t]d(Rg) (A) determined from the self-diffusion constant D, eq. ^, versus n for T = 448 K. 
Also shown for four values of n is the viscosity rj (•) measured in our non-equilibrium molecular 
dynamics simulations. 
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TABLE I. Intramolecular interaction parameters. 



Model A (UA a 



Model B (AUA b ) Model C (UA C 



Units 



bond length 
k& (bending) 

e b 

ao (torsion) 
ai 
a 2 
&3 
a i 
as 



1.54 
124.18 
114.0° 
2.007 
4.012 
0.271 
-6.290 



1.54 
124.28 
114.6° 
2.062 
4.821 
0.162 
-6.218 
-0.324 
-0.502 



1.54 
124.18 
114.0° 
1.736 
4.500 
0.764 
-7.000 



A 

kcal/ (mol rad 2 ) 
kcal/mol 



1 Intramolecular parameters for n-alkanes from Ref. [13|,14|. Torsional potentials are taken from 



Ref. |0|. 

intramolecular parameters from Ref. [16|; torsional potential (d) was used. 

c The torsional potential is the "modified" one from Ref. |ll|. The bending term used here is slightly 

different from the one used by Paul et al; see [19]. 

TABLE II. Lennard-Jones potential parameters. 



Model 



Group 



a (A) 



e (kcal/mol) 



d(A) 



A (UA a ) 
B (AUA b ) 
C (UA C ) 



CH 3 
CH 2 
CH 3 
CH 2 
CH 3 
CH 2 



3.930 
3.930 
3.527 
3.527 
4.010 
4.010 



0.227 
0.093 
0.238 
0.159 
0.227 
0.093 



0.275 
0.370 



a From Ref. g|g4|. 
b From Ref. fj||. 

c The value used in Ref. for sch 3 was 0.226 kcal/mol and the a for both CH3 and CH 2 united 
atoms was 4.009A 
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TABLE III. Summary of the constant density (p = 0.766 g/cm 3 ) simulation results for model 
A at T = 448 K. Results for the mean squared radius of gyration (R G ) and the end-to-end 
distance (R 2 ) are expressed in A 2 . The eigenvalues of the mass tensor (if) satisfy the equality 
(R%) = I 2 + I 2 + I 2 , with I 2 > l\ > I 2 . Results for the self-diffusion and rotational relaxation 
time tr are from equilibrium molecular dynamics simulations, while the results for the viscosity 
77 are from non-equilibrium molecular dynamics simulations. The total length of the run and the 
number of molecules used are also shown. Uncertainties in the last reported digit (s) are given in 
parenthesis. 



n 


N 


t (ns) 


(Rl) 


\\I(R 2 G ) \ 2 /(R G ) 


(R 2 ) 


D (10~ 6 cm 2 /sec) 


tr (ps) 


EMD 6 


64 


1.0 


4.39(1) 


0.897(1) 0.082(1) 


31.9(1) 


25.8(7) 


11 


10 


64 


1.5 


10.7(1) 


0.900(1) 0.082(1) 


87.5(2) 


22.7(7) 


33 


16 


64 


1.5 


23.5(1) 


0.876(1) 0.105(1) 


194(1) 


15.5(7) 


86 


24 


100 


1.0 


44.1(2) 


0.847(2) 0.130(2) 


353(3) 


11.2(5) 


176 


30 


100 


1.5 


61.8(4) 


0.833(3) 0.140(2) 


484(6) 


8.56(31) 


260 


36 


100 


2.0 


79.8(7) 


0.821(3) 0.148(3) 


605(10) 


6.40(34) 


360 


44 


100 


2.0 


104(1) 


0.804(4) 0.159(3) 


758(16) 


5.27(30) 


513 


66 


100 


14.0 


178(1) 


0.785(2) 0.170(2) 


1240(17) 


3.12(8) 


1140 


n 


N 


t (ns) 


T) (CP) 










NEMD 16 


100 


18.5 


0.967(70) 










30 


64 


18.0 


1.56(6) 










44 


100 


13.5 


2.30(5) 










66 


100 


10.1 


3.22(5) 
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TABLE IV. Summary of the simulations at the experimental density for P ~ 0.1 MPa and 
T = 400 K. For each chain length, data is presented from simulations employing models A, B, 
and C. Densities p are quoted in gm/cm 3 . The results for the mean squared radius of gyration 
(Rq) and the end-to-end distance (R 2 ) are expressed in A 2 . The eigenvalues of the mass tensor 
(I 2 ) satisfy the equality (Rq) = l 2 + l 2 + if, with I 2 > l\ > I 2 . The total length of the run and the 
number of molecules used are also shown. Uncertainties in the last reported digit (s) are given in 
parenthesis. 



Substance 


Model 


N 


t (ns) 


(R 2 e) 


il/{Rl) 


ll/{R 2 G ) 


(R 2 ) 


D 


D e xp 


n-CiQ 


A 


64 


2.0 


24.1(1) 


0.881(1) 


0.102(1) 


201(1) 


23.2(1.3) 


20(1) 


p = 0.6999 a 


B 


64 


2.0 


21.8(1) 


0.858(1) 


0.119(1) 


173(1) 


18.1(9) 






C 


64 


2.0 


23.3(1) 


0.876(1) 


0.105(1) 


192(1) 


18.5(1.0) 




n-C 2 4 


A 


100 


1.0 


45.3(3) 


0.852(3) 




367(5) 


12.9(9) 


8.23(8 


p = 0.7260 b 


B 


100 


1.0 


40.1(4) 


0.833(4) 




309(6) 


8.90(62) 






C 


100 


2.0 


43.3(2) 


0.845(2) 


0.131(2) 


344(3) 


9.11(28) 




n-C 30 


A 


64 


4.0 


64.0(4) 


0.838(2) 


0.137(2) 


507(5) 


8.62(47) 


5.2(3) 


p = 0.7421 c 


B 


64 


4.0 


55.3(4) 


0.820(3) 


0.148(2) 


415(5) 


5.30(30) 






C 


64 


4.0 


61.0(4) 


0.833(3) 


0.140(2) 


477(6) 


6.03(34) 




71-C44 


A 


100 


4.0 


111(1) 


0.815(3) 


0.152(3) 


833(15) 


4.52(21) 


2.4(2) 


p = 0.7570 d 


B 


100 


4.0 


92.3(1.2) 


0.794(4) 


0.165(3) 


654(15) 


2.49(12) 






C 


100 


4.0 


104(1) 


0.807(4) 


0.157(3) 


765(16) 


3.17(17) 





Experimental densities at ten temperatures were obtained from Ref. [3^,31,3^] and used to create 
a linear extrapolation to T = 400 K. 
b From Ref. @. 

c Obtained by a linear interpolation between experimental densities from Ref. p4| . 

d Calculated based on information from Ref. | {L8| . While no experimental reference is given in 

Ref. jl^], a comparison of the value used was made to the value obtained by a linear interpolation 



between data from Ref. The two values differed by < 1%. 
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